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ABSTRACT 

We analyse the stability and nonlinear dynamical evolution of power-law accretion disc mod- 
els. These have midplane densities that follow radial power-laws, and have either temperature 
or entropy distributions that are strict power-law functions of cylindrical radius, R. We employ 
two different hydrodynamic codes to perform high resolution 2D-axisymmetric and 3D simu- 
lations that examine the long-term evolution of the disc models as a function of the power-law 
indices of the temperature or entropy, the thermal relaxation time of the fluid, and the disc 
viscosity. We present an accompanying stability analysis of the problem, based on asymptotic 
methods, that we use to interpret the results of the simulations. We find that axisymmetric 
disc models whose temperature or entropy profiles cause the equilibrium angular velocity to 
vary with height are unstable to the growth of modes with wavenumber ratios \kR/kz\ » 1 
when the thermodynamic response of the fluid is isothermal, or the thermal evolution time 
is comparable to or shorter than the local dynamical time scale. These discs are subject to 
the Goldreich-Schubert-Fricke (GSF) or 'vertical shear' linear instability. Development of the 
instability involves excitation of vertical breathing and corrugation modes in the disc, with the 
corrugation modes in particular being a feature of the nonlinear saturated state. Instability is 
found to operate when the dimensionless disc kinematic viscosity v < 10"^, corresponding to 
Reynolds numbers Re = Hc^/v > 2500. In three dimensions the instability generates a quasi- 
turbulent flow, and the associated Reynolds stress produces a fluctuating efl'ective viscosity 
coefficient whose mean value reaches a ~ 6 x 10"^ by the end of the simulation. The evolu- 
tion and saturation of the vertical shear instability in astrophysical disc models which include 
realistic treatments of the thermal physics has yet to be examined. Should it occur on either 
global or local scales, however, our results suggest that it will have significant consequences 
for their internal dynamics, transport properties, and observational appearance. 

Key words: accretion discs - instabilities - methods: numerical,analytical - planetary sys- 
tems: protoplanetary disks 



1 INTRODUCTION 

Accretion discs play important roles in a broad range of astrophysi- 
cal phenomena. Protostellar discs orbiting young stars provide con- 
duits through which most of the mass accretes during the star for- 
mation process, and they are the nascent environments for plane- 
tary system formation. Mass transfer through Roche lobe overflow 
in close binary systems leads to highly energetic and time variable 
phenomena because of accretion through a disc onto compact ob- 
jects such as white dwarfs in cataclysmic variables, and neutron 
stars or black holes in low mass X-ray binaries. Quasars and active 
galactic nuclei are powered by disc accretion onto supermassive 
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black holes. Understanding the dynamics and evolution of astro- 
physical discs is key to understanding these and related phenom- 
ena. 



Since the early work of |Shakura & Sunyaev] ( |1973| ) and 
|Lynden-Bell & Pringle| ( [1974| , significant eff'orts have been made 
to understand the internal dynamics of discs, including mecha- 
nism(s) that allow them to transport angular momentum and ac- 
crete at their observed rates. Several ideas based on dynamical 
instabilities in non-magnetised flows have been explored, includ- 
ing the Papaloizou-Pringle instability for thick accretion tori driven 
by unstable non-axisymmetric wave modes (Papaloizou & Pringle 
T984t, convective instability ([Cameron & Pine||1973| |Lin & Pa- 



paloizoul 19801 [Ruden et al.|1988 
tational instability { [Toomre|1964[ 



Ryu & Goodman" 1992'), sravi- 



in & Pringle 1987, Papaloizou] 



|& Savonije|1991| ), the global and subcritical baroclinic instabilities 
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( |Klahr & Bodenhe imer'200?, To hnson & Gammie|2006[ [Petersen| 
|et al.|2007{|Lesur & Papaloizou 20 10|), and the vertical shear insta- 
bihty ( Urpin & Brandenburg 1998 , Urp in|2003| |Arlt & Urpin 2004 ) 
which is closely related to the Goldreich-Schubert-Fricke (GSF) in- 
stabiUty ( |Goldreich & Schubert|[T967[ |Fricke|[T968l ) developed in 
the context of differentially rotating stars. The presence of a weak 
magnetic field in the disc, however, leads to the development of 
magnetohydrodynamic turbulence driven by the magnetorotational 
instability (Balbus & Hawley|199T||Hawley~&Balbus 1991 ), and 
it is now generally accepted that this is likely to be the source of 
anomalous viscosity in most accretion discs in which the magnetic 
field is well-coupled to the gas. 



In this paper we present an extensive analysis of the hydrody- 
namic stability and nonlinear dynamics of disc models with power- 
law midplane density distributions, and either temperature or en- 
tropy profiles that are power-law function of R only, where R is the 
cylindrical radius. Our investigation employs high resolution 2D- 
axisymmetric and 3D hydrodynamic simulations and a linear sta- 
bility analysis based on asymptotic methods. Models with power- 
law temperature profiles, adopting locally isothermal equations of 
state, have been used extensively in the study of protostellar disc 
dynamics and disc planet interactions. These studies normally as- 
sume the disc is viscous (e.g . |Kley et al.||2001[ |Cresswell et aTT] 
|2007l[Fragner & Nelsonl2010','D'Angelo & Lubow'2010'), or mag- 
netised (e.g.[Fromang & Nelson 2006 , Beckwith et al. 201 1), but in 
this study we include neither viscosity nor magnetic fields. Adop- 
tion of radial variations in temperature or entropy in the models 
causes them to have angular velocity profiles that are a function of 
both radius and height, Q(R,Z). The height- variation of Q is often 
referred to as the thermal wind in studies of atmospheric dynamics. 



We find that disc models for which dQ/dZ + 0, and which ex- 
perience thermal relaxation on ~ dynamical time scales or shorter, 
are unstable to the growth of modes with \kRlkz\ » 1, where 
{]iR, kz) are the radial and vertical wavenumbers. This instabil- 
ity appears to be closely related to the GSF instability, as studied 
by [Urpin & Brandenburg] (1998] ), [Urpml ( [200^ and |Arlt & U^in| 
P004| ) in the context of accretion discs, and confirmed by our 
own linear stability analysis. Growth of the instability is favoured 
when the thermodynamic response of the gas is isothermal, of 
near-isothermal, although the strength of this dependence varies 
with model parameters with steeper thermal gradients display- 
ing a greater tendency toward instability. The one 3D simulation 
we presents suggests that nonlinear development of the instability 
leads to a turbulent flow whose associated Reynolds stress leads 
to an eff'ective alpha parameter a ~ 10~^, causing non-negligible 
outward angular momentum transport. 



This paper is organised as follows. In Sect. [2] we present the 
basic equations of the problem, and the disc models we examine. 
In Sect. |3] we discuss the hydrodynamic stability of rotating shear 
flows, and previous work in the literature relevant to the present 
study. In Sect. [4] we describe the numerical methods employed, and 
in Sect. [5] we present the results of the nonlinear simulations. A 
stability analysis of the problem using asymptotic methods is pre- 
sented in Sect. [6] and a reworking of the analysis of |Goldreich &| 
|Schubert| ( |1967| ) for a fully compressible fluid is given in the ap- 
pendix. We draw our conclusions and discuss ideas for future work 
in Sect. El 



2 BASIC EQUATIONS 

In this paper we make use of both cylindrical polar coordinates (R, 
0, Z) and spherical polar coordinates (r, 6, 0). We solve the conti- 
nuity, momentum and internal energy equations of hydrodynamics 



dtipY) + V • [pVV] 

dt{e) + V • (^v) 



0, 

-VP - VO, 
-PV x + S-Q 



(1) 



where p is the density, v is the velocity, P is is the pressure, e is 
the internal energy per unit volume, S and Q are energy source and 
sink terms, and O = -GM/r is the gravitational potential due to the 
central star. Here G is the gravitational constant and M is the mass 
of the star. 



2.1 Disc models 

We are concerned primarily with two basic equilibrium disc models 
in this paper. In the first, the temperature, T, and midplane density, 
Pmid, are simple power-law functions of cylindrical radius: 



T(R) 



PmidW 



(2) 



(3) 



where Tq is the temperature at the fiducial radius Rq, and po is the 
midplane gas density there. Adopting an ideal gas equation of state 



-Tp, 



(4) 



where 9^ is the gas constant and is the mean molecular weight, we 
note that the isothermal sound speed is related to the temperature 
through the expression = jfi, such that q also represents the 
radial power-law dependence of cJ(P): 



cliR) ■■ 



(5) 



In the second disc model, we adopt a power-law function for 
the midplane density, as in eqn. (|3j, and specify the entropy func- 
tion, Ks, as a strict power-law function of R in the initial model: 



where the entropy function is defined through the expression 
P = Ky, 



(6) 



(7) 



and y is assumed to be constant. The entropy per unit mass is given 
by 



5 =Cviog(rp'-'') 



(8) 



where Cy is the specific heat at constant volume, and the entropy 
function is given in terms of the entropy by the expression 



Ks = Cy(y - l)exp — 



(9) 



2.1.1 Equilibrium solutions 

In order to construct initial conditions for our simulations we need 
to obtain equilibrium disc models. The equations of force balance 
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in the radial and vertical directions are given by 
GMR 



(R^ + Z2)3/2 
GMZ 



pdR 



1 dP 



= 0. 



(10) 



(11) 



(R^ + Z2)3/2 p dZ 

Combining eqns. (To| , ([TTJ, (|2]) and ([3]) leads to expressions for the 
equilibrium density and angular velocity, Q, as functions of (R, Z) 
for the disc with a power-law temperature profile: 

IGM 



p{R,Z) 



Po 



1 



a(R,z) 



4l^ 



V^^Tz2 
+ (1 + ^)- 



(12) 



qR 



V^^Tz2 



1/2 



(13) 



where Qk = ^|GM|R^ is the keplerian angular velocity at radius 
R, and // = cJCIk is the local disc scale height (also see |Takeuchi| 
|& Liri[2002l [Fromang et al.|2011a| ). The definition of given in 
eqn. ([5]) implies that 

/ n \(^+3)/2 

H-H.,(-) (14) 



where Hq - cqI ^GM/R^ is the disc scale height at radius Rq. 

Similarly, the equilibrium density and angular velocity for the 
disc model with a power-law entropy function profile are given by: 

^1/(7-1) 

(15) 



p(R,Z) 



Q(R,Z) 



{y-\)GM 



1 



V/?2 + Z2 



+ (1 + 5) - 



5/? 



1/2 



(16) 



yM^ ■ " ■ " V/?2 + Z2 

where Almid = ^Klcirmd is the Mach number at the disc mid- 
plane and M - vk/ as is the Mach number at each disc location. 
vk = yjGM/R is the keplerian velocity at radius R, and the adia- 
batic sound speed = yJyPjp, which takes the value amid at the 
disc midplane. We note that combining equations (|4]), ([6|, ([7| and 
|T5] l demonstrates that the temperature in this model is a function 
of both R and Z in general: 



T{R,Z)^K,p 



n' 



(17) 



As such, this model provides a useful contrast to the one with tem- 
perature constant on cylinders, and is convenient to implement nu- 
merically because of the existence of analytic solutions for the equi- 
librium p{R, Z) and Q(/?, Z) profiles. 

Equations (|2|, ([3]), |T2] ) and ([13]) fully specify the initial disc 
models with power-law temperature profiles that we examine in this 
paper, subject to appropriate choices for p and q. The expressions 
([3]) and ([6|, along with ( p3] l and ([T6|, specify the initial models for 
which the entropy function is a power-law function of cylindrical 
radius, subject again to appropriate choices for the power-law ex- 
ponents p and s. We note that for all models in which the initial 
temperature, T, or entropy function, Ks, are strict power-law func- 
tions of R, the equilibirum angular velocities are explicit functions 
of both R and Z, a fact that appears to play a key role in the disc 
instability that we examine in this paper. The dependence of Q on 
Z is often referred to as the 'thermal wind' in studies of planetary 
atmosphere dynamics. 



2.2 Thermodynamic evolution 

The thermodynamic evolution of both disc models described above 
in Sect. |2.1| is assumed to be one of three types in the simulations 



presented here: locally isothermal, for which the local temperature 
at each {R, Z) position in the disc is kept strictly equal to its orig- 
inal value; isentropic, where the entropy of the fluid is kept con- 
stant (equivalent to there being no source/sink term in the energy 
eqn. |[T)); thermally relaxing, where we relax the temperature at 
each location in the disc toward its initial value on some time scale, 
Treiax- The thermal relaxation model we adopt is 



dT 
d7 



{T - To) 

'^relax 



(18) 



where Tq is the initial temperature. For simplicity, we assume that 
Treiax is a function of R, being a fixed multiple or fraction of the 
local keplerian orbital period. Equation (T8| ) has a simple analytic 
solution of the form 



T{t + AO = To + {T{t) - To) exp 



A^ 

'^relax / 



(19) 



where T(t) is the temperature at time t, and T(t + AO is the temper- 
ature at some later time t + A^. 

In the locally isothermal models we use an isothermal equa- 
tion of state P - c^p, and do not evolve the energy equation in 
eqn. ([TJ. In the isentropic models we use the equation of state 
P - {y - \)e, solve the energy equation in eqn. ([TJ, and neglect 
the source and sink terms. The energy equation is also solved in the 
thermally relaxing models, along with eqn. ( [Tsj l which plays the 
role of the source and sink terms in the energy equation ([TJ. 



3 HYDRODYNAMIC STABILITY OF DISC MODELS 

3.1 The Rayleigh and Solberg-H0iland criteria 

The Rayleigh criterion indicates that accretion discs with strictly 
keplerian angular velocity profiles are hydrodynamically stable 
since 



dl 
dr 



>0, 



where j - R^Q and Q = ^]GM|R^ . More generally, a diff'er- 
entially rotating, compressible fluid with angular velocity varying 
with height and radius, Q(/?,Z), subject to axisymmetric isentropic 
perturbations (i.e. DS /Dt = 0, where D/Dt is the total time deriva- 
tive for fluid elements) is stable if both of the Solberg-H0iland cri- 
teria are satisfied (e.g. |Tassoul|1978j ): 

1 df 



dP_(dj^ds_ 

dZ\~dRdZ 



~dZ~dR 



>0. 



(20) 



(21) 



Accretion disc models generally possess negative radial pressure 
gradients, and a negative vertical pressure gradient in the disc hemi- 
sphere above the midplane, leading to stability criteria ( [20j and (|2TJ 
in the form: 



1^ 
R^ OR 



1 

pCp 

dj^ds_ 

~dR~dZ 



( dP 


ds 


dP 


ds\ 


\ OR 


dR^ 


dz 


dz) 



>0 



'dZdR 



>0. 



(22) 



(23) 



Considering a disc with a strictly keplerian j profile we see that 
stability according to eqn. ([23j requires dS /dZ > 0, in agreement 
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with the Schwarzschild condition for convective stabiUty. Condi- 
tion ( [22| shows that a large ampUtude negative radial entropy gra- 
dient dS IdR < can also drive instability in principle, provided the 
gradient is strong enough to overcome the positive angular momen- 
tum gradient. This is distinct from the global baroclinic and/or sub- 
critical baroclinic instability discussed in the introduction that does 
not require violation of the Solberg-H0iland criteria to operate, but 
does require thermal evolution of the fluid on time scales similar 
to the dynamical time to re-establish the initial radial entropy gra- 
dient, in addition to the presence of finite amplitude perturbations 
( |Lesur & Papaloizou|20 10^. 

Considering a quasi-keplerian disc whose angular velocity de- 
pends on height and radius, 0.{R, Z), eqn. ( [23] ) shows that such a 
disc with dS IdZ > can be destabilised through the combination 
{dfldZ){dSldR)>0. 

We consider two basic disc models in this paper that are the 
subject of nonlinear simulations in which the fluid evolution is isen- 
tropic and for which the Solberg-H0iland criteria determine hy- 
drodynamic stability. One model assumes the temperature profile 
is a strict power law function of radius, R, such that the density 
and angular velocity are given by eqns. |T2| and |T3] ). With val- 
ues p - -1.5 and q - -\, -0.5, -0.25 and in eqns. ([2| and ([3]), 
respectively, these discs are stable according to criterion ( [22| . The 
term involving dS IdR provides the only destabilising contribution 
(because dS jdZ > 0), but is always far too small to overcome the 
positive radial gradient in f. This disc is also stable according to 
criterion ([23} as (df/dR)(dS/dZ) > (df/dZ < 0)(dS/dR). 

The second set of disc models we consider assume that the 
entropy is a strict power law function of radius R. The density and 
angular velocity are given by the expressions ( p3] l and ([16}, and 
the values (p = 0, s = -\) and (p - -1.5, s - 0) are adopted 
in eqns. ([3}and ([6}. Both of these models are stable according to 
criterion ( [22| ) as the destabilising term proportional to dS I OR is 
either zero or too small (note that dS jdZ - in these models). 
Criterion ([23}, however, predicts instability for the (/? = 0, ^ = -1) 
model because only the second term on the left-hand side is non 
zero, with df jdZ < and dS jdR < 0. The isentropic simulation 
that we present in the results section later in this paper will provide 
an example of nonlinear evolution of this model that is unstable 
according to one of the Solberg-H0iland criteria. The model with 
{p - -\.5, s - Qi) is stable according to criterion ([23}. 

3.2 The Goldreich-Schubert-Fricke instability 

When thermal and viscous diff'usion play a role the stability of ro- 
tating flows are no longer controlled by the Solberg-H0iland crite- 
ria, but instead are determined by stability criteria obtained origi- 
nally by Goldreich & Schubertl(|1967|) and Fricke ( 1968 ) in applica- 
tion to the radiative zones of difl'erentially rotating stars. Axisym- 
metric rotating flows are susceptible to the Goldreich-Schubert- 
Fricke (GSF) instability when viscous difl'usion is much less ef- 
ficient than thermal diff'usion, such that a fluid element retains its 
initial angular momentum but quickly attains the entropy of the sur- 
rounding fluid when perturbed from its equilibrium location. Under 
these circumstances the stabilising influence of entropy gradients 
provided by the Solberg-H0iland criteria diminish and instability 
ensues for wave modes satisfying the instability criterion 

|^-^|^<0. (24) 
dR kzdZ ^ ^ 

For a rotating flow where the angular momentum is a function of 
both R and Z, and the appropriate conditions on thermal and viscous 



difl'usion are satisfied, unstable modes are guaranteed to exist since 
wavevectors with ratios kRikz that satisfy eqn. ([24} can always be 
found. In general, we expect \d j/dR\ » \d j/dZ\ in a quasi-keplerian 
accretion disc, such that unstable modes will have \kR/kz\ » 1 (i.e. 
unstable modes will have radial wavelengths that are very much 
shorter than vertical ones). 

Application of the GSF instability to accretion discs has not 
received a great deal of attention in the literature (but see the dis- 
cussion below). Indeed, the study presented in this paper has arisen 
from unrelated attempts to generate 3D models of protoplanetary 
disc models within which turbulence is generated by the MRI in 
active layers near the disc surface, and in which extensive dead 
zones that are magnetically inactive (or at least stable against the 
MRI) exist near the midplane. The adoption of a locally isothermal 
equation of state in these models, combined with the absence of 
a physical viscosity, rendered them unstable to the growth of ver- 
tical corrugation oscillations that in the nonlinear regime became 
quite violent. This behaviour appears to be a nonlinear manifesta- 
tion of the GSF instability, which we study here in more detail. The 
two classes of models we consider (temperature constant on cylin- 
ders and entropy constant on cylinders) both have angular velocity 
profiles that depend on both R and Z, and we expect that models 
in which perturbations evolve quasi-isothermally will display the 
GSF instability. 

Previous analysis of the GSF instability in the context of 
accretion discs was initiated by jUrpin & Brandenburg (1998 ) in 
which they presented a local linear analysis utilising the short- 
wavelength approximation. A more extensive analysis was pre- 
sented in |Urpin| ( [2003 ) where it was suggested that the GSF may 
provide a source of hydrodynamic turbulence in accretion discs, 
and nonlinear simulations were presented by Arlt & Urpin ( 2004|). 
These nonlinear simulations adopted a basic disc model with a 
strictly isothermal equation of state. As such, the underlying equi- 
librium disc model has an angular velocity profile that is indepen- 
dent of height, Z, (seeeqnJT3} and should not be susceptible to the 
GSF instability. Arlt & Urpin| ( [2004] ), however, adopted initial con- 
ditions that allowed relaxation around the equilibrium state leading 
to variations of Q with height, and when initial perturbations with 
\kR/kz\ » 1 were applied growth of the GSF was observed. It is 
also worth noting the related linear and nonlinear study of isen- 
tropic accretion disc models presented by [Riidiger et al.| ( |2002| ) at 
this point. This study examined the stability of discs whose angu- 
lar momentum and entropy profiles rendered them stable according 
to the previously discussed Solberg-H0iland criteria. As expected, 
applied perturbations were found to always decay. We present a 
linear analysis of the instability we study in Sect.[6]for a fully com- 
pressible fluid under the assumption that the equation of state is 
locally isothermal, and in the appendix we present an analysis of 
the problem that closely follows the derivation of the GSF stability 
in .Goldreich & Schubert] ( [T967| ). 



4 NUMERICAL METHODS 

The simulations presented in this paper were performed using two 
difl'erent codes that utilise very difl'erent numerical schemes. We 
use an older version of nirvana, which uses an algorithm very sim- 
ilar to the ZEUS code to solve the equations of ideal MHD ( Zieglerj 
|& Yorkel 19971 [Stone & Norman] 1992| ). This scheme uses operator 
splitting, dividing the governing equations into source and transport 
terms. Advection is performed using the second-order monotonic 
transport scheme ( [van Leer|1977| ). We also use the more modem 
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Table 1. Simulation parameters and results. Labels beginning with a 'T' 
denote runs with T being a function of R. Those with a 'K' denotes runs 
where Kg = Ks(R). The letter R denotes reflecting boundary conditions, 
and a letter O denotes outflow b.c's. Digits after the hyphens denote thermal 
relaxation times. Note that all runs have Hq/Rq = 0.05. 



Simulation 



J or S TRelax 



Unstable 



7 


TlR/O-0 


-1.5 


-1 


0.00 


1328 X 1000 X 1 


Y 


T2R/O-0 


-1.5 


-0.5 


0.00 


1328 X 1000 X 1 


Y 


T3R/O-0 


-1.5 


-0.25 


0.00 


1328 X 1000 X 1 


Y 


T4R/O-0 


-1.5 





0.00 


1328 X 1000 X 1 


N 


T5R-0.01 


-1.5 


-1 


0.01 


1328 X 1000 X 1 


Y 


T6R-0.1 


-1.5 


-1 


0.10 


1328 X 1000 X 1 


N 


T7R-1.0 


-1.5 


-1 


1.00 


1328 X 1000 X 1 


N 


T8R-10.0 


-1.5 


-1 


10.0 


1328 X 1000 X 1 


N 


T9R-00 


-1.5 


-1 


oo 


1328 X 1000 X 1 


N 


KlR-0 





-1 


0.00 


1328 X 1000 X 1 


Y 


K5R-0.01 





-1 


0.01 


1328 X 1000 X 1 


Y 


K6R-0.1 





-1 


0.10 


1328 X 1000 X 1 


Y 


K7R-1.0 





-1 


1.00 


1328 X 1000 X 1 


Y 


K8R-10.0 





-1 


10.0 


1328 X 1000 X 1 


Y 


K9R-CO 





-1 


oo 


1328 X 1000 X 1 


N 


KlOR-0.01 


-1.5 





0.01 


1328 X 1000 X 1 


N 


T1R-0-3D 


-1.5 


-1 


0.00 


1328 X 1000x300 


Y 



NiRVANA-iii code, which is a second order Godunov-type MHD code 
( |Ziegler|[2004| ), which has recently been extended to orthogonal- 
curvilinear coordinate systems (Ziegler 2011 ). All presented simu- 
lations were performed using a standard spherical coordinate sys- 
tem (r, 6, 0). 

4.1 Initial and boundary conditions 

The numerical study presented here applies to the two general 
classes of disc model discussed in Sect. |2.1| and as such we adopt 
a numerical set up that is not specific to any particular physical 
system (although our motivation for undertaking this study arose 
from earlier attempts to establish stationary equilibrium solutions 
for protostellar disc models). Based on numerous test calculations 
performed with a wide range of resolutions during an early stage of 
this project, we know that the instability we study here is charac- 
terised by having a radial wavelength much shorter than the vertical 
wavelength (i.e. kR » kz) during its early growth phase. Conse- 
quently we consider disc models with fairly narrow radial domains 
to facilitate high resolution simulations. The spherical polar grid 
we adopt has inner radius = Rq = \ and outer radius rout = 2, 
and most simulations we present are axisymmetric. The one non- 
axisymmetric simulation we present covers a restricted azimuthal 
domain of 7t/4, again to facilitate a high resolution study. 

Disc models in which the initial temperature is a strict function 
of cylindrical radius, R, have meridional domains 7t/2 - SHq/Rq < 
6 < 7t/2+5Ho/Ro- For a disc with radial temperature profile q = -\, 
corresponding to a disc with constant H/R, the meridional domain 
covers +5 scale heights above and below the midplane. For larger 
(less negative) values of q, the disc covers +5 scale heights at the 
inner radius, but a reduced number of scale heights as one moves 
out in radius. Prior to initiating these simulations, the disc models 
are specified using eqns. (jSj, ([2|, (T2| and (\3\ . In all models the 
initial velocity field was seeded with random noise distributed uni- 



formly in each component with a peak amplitude equal to 1 % of 
the local sound speed. 

Discs for which the initial entropy function, K^, depends only 
on R have physical surfaces where the density goes to zero. In these 
models the meridional boundaries are placed at a location that is 
5% smaller than the angular distance from the midplane where the 
density vanishes. The initial models are specified using eqns. (51, 
([6|, (T5| l and ([16]). In order to determine the value of in eqn. ( 6J, 
we specify the midplane Mach number, M^id^ at radius Rq. We 
adopt Mjnid = 20 to be consistent with the models for which HjR- 
0.05. Seed noise with amplitude equal to 1 % of the sound speed was 
again added to all velocity components. 

For most simulations we adopt standard outflow or reflecting 
boundary conditions at the inner and outer radial boundaries. Pe- 
riodic boundaries are applied in the azimuthal direction, and ei- 
ther standard outflow or reflecting conditions are applied at the 
meridional boundaries (all simulations performed using nirvana- 
III adopted outflow boundary conditions at the radial and merdional 
boundaries). The density is obtained in ghost zones by means of lin- 
ear extrapolation. A variety of diff'erent boundary conditions were 
used in test simulations at an early stage of this project, and the 
results were found to be insensitive to the choice adopted. These 
tests included the adoption of damping boundary conditions that 
absorb incoming waves, indicating that reflecting boundaries are 
not required to drive the instability discussed in this paper. 

Owing to the unsplit character of finite volume schemes, such 
as used in nirvana-iii, it is difficult numerically to preserve static 
equilibria ( [Zingale et aL|2002| . In particular, problems arise with 
the constant extrapolation of the density profile in the vertical di- 
rection. This is because the weight of the gas is not balanced by the 
now vanishing pressure gradient in the case of an isothermal equa- 
tion of state, leading to a standing accretion shock in the first grid 
cell of the meridional domain. In the case of an adiabatic evolution 
equation, the same problem arises in the boundary condition for the 
thermal energy (which enters the boundary condition for the total 
energy). To alleviate this situation, we obtain the density (and in 
the case of solving an energy equation, the pressure) by integrat- 
ing the hydrostatic equilibrium in each cell adjacent to the domain 
boundary. This is done with a second-order Runge-Kutta shooting 
method. 

Unless indicated otherwise, we use a fixed resolution of A^^ x 
of 1328 X 1000 grid cells for our nirvana runs with Hq/Rq - 0.05 
or Almid = 20 at r = i?o- The nirvana-iii runs used a resolution of 
1344 X 1024 

We adopt a system of units in which M - \,G - 1 and Rq - \. 
When presenting our results, the unit of time is the orbital period at 
the disc inner edge, = In. 



5 RESULTS 

The main aims of the following simulations are to delineate condi- 
tions under which the disc models outlined in Sect. 12. II are unsta- 
ble to the growth of the vertical shear instability, and to examine 
the effect that different physical and numerical set-ups have on the 
growth and evolution of this instability. We also aim to characterise 
the final saturated state of these unstable discs, although our adop- 
tion of mainly axisymmetric simulations provides some restriction 
in achieving this final aim. 

We define the volume integrated meridional and radial kinetic 
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Figure 1. Time evolution of the normalised perturbed kinetic energy in the 
meridional and radial coordinate directions for model TlR-0 with p = -1.5, 
q = -I and reflecting boundary conditions at the meridional boundaries. 

energies through the expressions 

ee^^j^pvldV, er^^j^pv^dV. (25) 

When presenting our results we normalise these energies by the to- 
tal kinetic energy contained in keplerian motion in the disc initially. 

We begin discussion of our results below by describing sim- 
ulations of locally isothermal discs for which T{R) ~ and 
p(R) ~ R~P. We describe one fiducial model in detail, before dis- 
cussing briefly the influence of the temperature profile in control- 
ling the instability. We present a comparision between results ob- 
tained using the two codes described in Sect. |4] and also demon- 
strate how the instability evolves as a function of disc viscosity. 

The next set of results we present are for disc models in which 
the initial temperature is a strict function of R, but we set y = 1 .4, 
solve the energy in eqn. ([T]), and allow the temperature to relax 
toward its initial value using eqn. |T8] l. In this section, we exam- 
ine how the thermal relaxation rate contols the instability, covering 
the full range of thermodynamic behaviour from locally isothermal 
through to isentropic. 

In the penultimate part of our numerical study, we consider 
disc models for which the entropy function, is a strict function 
of R, and again employ thermal relaxation to examine the condi- 
tions under which accretion discs display the vertical shear instabil- 
ity. The final numerical experiment we present examines the insta- 
bility in a full 3D model, and provides an estimate of the Reynolds 
stress induced by the instability. 

We present an analytic model in the discussion section which 
illustrates the basic mechanism of the instability and delineates the 
conditions under which it operates. 

5.1 A fiducial model 

We begin presentation of the simulation results by discussing one 
particular model in detail to illustrate the nature of the instability 
that is the focus of this paper. The fiducial model is TRl-0 listed in 
Table [T] with temperature, T, defined as a function of R only, and 
p = -1.5 and q = -1 in eqns. ^ and ([2]). A locally isothermal 
equation of state is adopted. As such, this disc model has param- 
eters very similar to those used in numerous previous studies of 
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Figure 4. Time evolution of the vertical centre of mass position for each 
radial location in the disc for the fiducial model TlR-0 with H/R = 0.05. 
Note that the vertical centre of mass position has been normalised by the 
local scale height at each radius. Starting from bottom to top the plots cor- 
respond to times (in orbits): 9 x 10"^, 9.18, 27.86, 65.00, 92.66. The multi- 
plicative factor indicated in each legend causes the maximum amplitude of 
the normalised c.o.m. position in each graph to equal unity. 

disc related phenomena (e.g. |Kley et al. [2001 [[Cress well & Nelsonj 
[2006l[F7omang et al.|20l7bl[Pierens & Nelson|2010[ ), although we 
focus primarily on the inviscid non-magnetised evolution here. 

The time evolution of the normalised meridional and radial 
kinetic energies defined in eqn. ( [25] ) are shown in Fig. [T] The in- 
tial values at ^ = originate from the seed noise, and we observe 
that after ~ 10 orbits, during which the perturbed kinetic energies 
damp slightly, rapid growth of the perturbation energies arises. The 
normalised energies reach non linear saturation after ~ 400 orbits 
having reached values of a few xlO~^. Inspection of the evolution 
of the sum of the meridional plus radial kinetic on a log-linear plot 
indicates that the linear growth rate of the perturbed energy in Fig.[T] 
is ^ 0.24 orbit"^ 

Contour plots of vertical velocity perturbations, Vz, that arise 
at diff'erent stages of the disc evolution are shown in Figs. |2] and 
|3] These two figures show the perturbed velocity field at identical 
times, but whereas Fig.[2|maps linearly between the velocity values 
and the grey-scale, Fig^plots the values sign(vz) x Ivzl^^"^ so that 
the grey- scale is stretched to enable the morphology of the pertur- 
bations to be more clearly discerned. Both figures demonstrate that 
perturbations start to grow near the upper and lower disc surfaces, 
where |dQ/dZ| is largest, and toward the inner edge of the disc. 
The perturbations are characterised by having short radial and long 
vertical wavelengths, as expected for the vertical shear instability 
described in Sects. [3.2[ and [6] The short radial wavelength gives 
rise to significant radial shear in the vertical velocity dvz/dR, and 
this apparently causes small scale eddies to form at the shearing 
interfaces. As time proceeds the instability extends toward the disc 
midplane and out to larger radii, until the entire disc participates in 
the instability (although it should be noted that the midplane where 
dQ/dZ = is formally stable to local growth of the vertical shear 
instability). 

Close inspection of the lower panels in Fig. [3] show that the 
velocity perturbations have odd symmetry about the midplane ini- 
tially (this is particularly apparent in the fourth panel between radii 
1.1 < r < 1.4). In other words the disc exhibits a breathing mo- 
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Figure 2. Edge-on contours of the perturbed vertical velocity as a function of R, Z and time for model TlR-0. 
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Figure 3. Edge-on contours of the perturbed vertical velocity as a function of R, Z and time for model TRl-0. Note that for clarity, the grey-scale of the image 
has been streteched by plotting the quantity sign(vz) x Ivzl^^"^. Note that the spectrum bar shows values of v^J^. 



tion about the midplane as the instability first becomes apparent at 
lower disc latitudes. As time proceeds, however, the velocity pertur- 
bations become symmetric about the midplane as demonstrated by 
the fifth and sixth panels. These perturbations correspond to cor- 
rugation of the disc characterised by coherent oscillations of the 
vertical centre of mass position whose phase depends on radius in 
a time dependent manner. The development of this disc corrugation 
is illustrated by Fig. |4| which shows the vertical centre of mass po- 



sition of the disc at each radius for five diff'erent times (note that 
each plot is off"- set in the vertical direction to aid clarity, and each 
curve has been multiplied by a unique factor so that the corruga- 
tion may be observed). The vertical centre of mass has been nor- 
malised by the local disc scale height at each radius. Moving from 
the lower curve to the upper curve, we note that the vertical centre 
of mass position has a very small variation with radius after one 
time step, but this becomes progressively larger in amplitude and 
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Figure 5. Edge-on contours of the disc relative density perturbations 6p/po as a function of R, Z and time. Note that we have effectively stretched the grey-scale 
by plotting the quantity sign((5p) x ^jSplpQ 



more spatially coherent as time progresses (times corresponding to 
each curve are given in the figure caption). The final curve, cor- 
responding to an evolution time of ~ 92 orbits, has a maximum 
vertical displacement approximately equal to 0.006//. It is interest- 
ing to note that the initial disc instability begins with \kRlkz\ » 1 
due to the short radial wavelengths of the fastest growing modes 
of the vertical shear instability. As the disc approaches the nonlin- 
ear state, however, the development of coherent corrugation waves 
causes the radial wavelengths of the most apparent perturbations to 
approach or modestly exceed the local scale height. At the end of 
the simulation (~ 420 orbits) the maximum vertical displacement 
of the disc centre of mass reaches ~ 0.01//, but we are cautious 
about interpreting the results at this late stage of evolution as the 
reflecting boundary conditions may play a role. 

Contour plots of the density perturbations 6plpQ correspond- 
ing to the previously discussed velocity contours are displayed in 
Fig.|5] As with the velocity contours, we see perturbations first arise 
near the disc surface. Coherent density structures first become vis- 
ible in the third panel, but are substantially more apparent in the 
lower panels. 

The final saturated state consists of locally unstable disc an- 
nuli that oscillate vertically at close to the local frequency, 0.\{R), 
superposed on which are a spectrum of oscillations with diff'erent 
frequencies caused by travelling waves excited by vertical oscilla- 
tions at neighbouring disc radii. A region of the disc lying at inter- 
mediate radii will thus experience a locally generated corrugation, 
in addition to inward travelling corrugation waves that propagate 
as inertial modes (or r-modes) launched from exterior disc loca- 
tions, and outward propagating corrugation waves propagating as 
acoustic or fundamental modes launched from interior disc loca- 
tions ( |Lubow & Pringle|1993|[Lubow & Ogilvie|1998{ ). 

Although we only present simulations with finite amplitude 
initial perturbations to the velocity fields in this paper, we have 
conducted numerous experiments in which the peak amplitude of 



the imposed perturbations varies, including cases where perturba- 
tions just grow from numerical round-off" errors. Although this re- 
quires a larger time interval for the instability to become apparent, 
we nonetheless observe instability for all perturbation amplitudes, 
demonstrating that the instability is linear. 

5.2 Evolution as a function of the radial temperature profile 

In this section we discuss simulations TlR-0 to T4R-0 (which 
utilise reflecting boundary conditions in the meridional direction) 
and TlO-0 to T4O-0 (which use open boundary conditions). These 
simulations have p - -1.5 and a range of values for q running from 
^ = -1 (a constant HjR disc) up to ^ = (a purely isothermal disc 
in which HjR oc R^'^). 

The left panel of Fig. [6] shows the time evolution of the nor- 
malised perturbed kinetic energy summed over the radial and mer- 
dional directions in the disc models TlR-0 to T4R-0. We see that 
as the value of q increases from q - -I through to ^ = -0.25, the 
growth rate of the instability decreases, and for ^ = it switches 
ofl" altogether. Inspection of the evolution of the sum of the merid- 
ional and radial kinetic energies on a log-linear plot indicates that 
the linear growth rate for the q - -0.5 case is ^ 0.12 orbit~\ and 
for the q - -0.25 run is ^ 0.052 orbit"^ (to be contrasted with 
the growth rate ^0.24 orbit" ^ obtained for the ^ = -1 run). The 
azimuthal velocity in the ^ = model is independent of Z, as in- 
dicated by eqn. ([13]), so the observed stability of this disc model is 
in agreement with the expectations discussed in Sect. |3] It is also 
noteworthy that the saturated values of the perturbation energies in 
each unstable disc, normalised to the total energy in keplerian mo- 
tion, difl'er substantially from one another in accord with the trend 
in the growth rates. 

The right panel of Fig. [6] shows the corrugation energies for 
runs TlO-0 to T4O-0, which are identical to runs TlR-0 to T4R- 
except that the boundary conditions applied at the meridional 
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Figure 6. Time evolution of the perturbed meridional plus radial ki- 
netic energies (normalised) as a function of the radial temperature profile. 
The left/right panel shows results for simulations which adopted reflect- 
ing/outflow boundary conditions at the upper and lower disc surfaces. 



boundaries are outflow rather than reflecting. During the growth 
phase of the instability the results of the TlO-0 - T4O-0 runs are 
essentially identical to the corresponding TlR-0 - T4R-0 calcula- 
tions, and the saturated energies are very similar. We also observe 
the important result that the transition between stable and unstable 
disc models is independent of the boundary conditions, with both 
q = models showing decay of the perturbed meridional and radial 
kinetic energies during the simulations. It is clear that the existence 
of a radial temperature profile plays a fundamental role in deter- 
mining whether or not a disc becomes unstable. 



5.3 Code comparison 

Numerical modelling of the Navier-Stokes equations can pose 
arcane pitfalls, particularly if hydrodynamic instabilities are in- 
volved. In the absence of rigorous analytical reference solutions, 
it has become customary to substantiate the physical reliability of 
the solutions obtained by means of code comparisons. While we 
have already shown the generality of a GSF-like instability occur- 
ring under various physical settings, we demonstrate here its com- 
parative development using two numerical schemes. Although very 
similar in their names, the two codes we used are fundamentally 
different in the numerical schemes they apply, nirvana utilises the 
same finite-difference scheme as the zeus code, whereas nirvana-iii 
applies a finite-volume Godunov scheme very similar to the ones 
used in, e.g., ramses ( |Teyssier|2002l ) or ATHENA ( [Gardiner & Stone] 
[2008] ). 

We have run the simulations labelled as TlO-0 to T4O-0 in 
Table 1 using both codes. The codes generally agree well in the 
development of the instability and its saturation level as demon- 
strated by Fig. [7] although it should be noted that the realisation of 
the initial seed noise in the two runs is slightly different. The de- 
cay of these initial seed perturbations occurs slightly faster in the 
NIRVANA-III simulations, and the saturation amplitude is slightly 
smaller. A detailed look at the vertical velocity perturbations at time 
t ^ 12.8 orbits is shown in Fig. [8] and both codes show the char- 
acteristic \kR/kz\ » 1 perturbations associated with the early de- 
velopment of the vertical shear instability. The codes are in decent 
agreement about the magnitude of the velocity perturbations and 
also in the dominant wavelength of instability. The codes, however, 
also show some differences in their solution at this time. The con- 
tinued presence of the initial seed noise is more apparent in the 
NIRVANA run than in the nirvana-iii run, in agreement with the evo- 
lution of the kinetic energies in Fig. |7] and the nirvana run shows 
a greater degree of structure in the velocity perturbations, perhaps 
indicative of higher-order modes being present at this time. Overall 
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Figure 7. Comparison of the time evolution of the normalised (meridional 
+ radial) kinetic energies for models TlO-0 to T4O-0 that were run with 
the two codes. 
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Figure 8. Comparison between the vertical velocity perturbations at time 
t = 12.8 orbits for the nirvana and nirvana-iii runs. The results are for sim- 
ulation TlO-0. 



the comparision is very satisfactory and demonstrates that the con- 
ditions for the vertical shear instability to occur are predicted ac- 
curately by both numerical schemes, which also show reasonable 
agreement for the growth rates under different radial temperature 
profiles. 



5.4 Evolution as a function of viscosity 

In this section we present results from simulations that examine the 
amplitude of the saturated state as a function of imposed viscos- 
ity. We apply a constant kinematic viscosity, y, to the disc model 
TlR-0 and vary its value between 10~^ < y < 10~^ (a value of 
y = 10~^ corresponds to the S hakura-Sunyaev viscous st ress pa- 
1 ( |Shakura & Sunyaev|l973] ) 



rameter a = 4x10 BiR 



and to 



a Reynolds number Re = Hcjv - 2500). The results are shown in 
Fig. [9] which shows the time evolution of the perturbed meridional 
plus radial kinetic energies. As expected, the results have a strong 
dependence on viscosity. For y = 10~^ the instability is damped 
completely, which explains why previous 3D simulations of locally 
isothermal discs have not reported seeing the vertical shear insta- 
bility (e.g.|Kley et al.|[200T| [Cress well & Nels6n||2006l |Fromang| 



|et al.|2011b[|Pierens & Nelson|2010| ). For decreasing values of 1 
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Figure 9. Time evolution of the normalised perturbed kinetic energy in 
the meridional and radial coordinate directions for model TlR-0 with 
p = -1.5, q = -1 and reflecting boundary conditions at the meridional 
boundaries. Each curve corresponds to a difl'erent value of the imposed 
kinematic viscosity, v, as indicated. 

the amplitude of instability increases, until at a value of y = 10~^ 
there is little difference between the result in Fig.|9]and the inviscid 
result shown in the left panel of Fig [6] 

Interestingly, it is found that in fully turbulent models where 
the MRI is active throughout the disc and a ^ 0.01, the corruga- 
tion instability is not observed ( [Fromang & Nelson|2006|). We have 
computed models similar to those presented in jFromang & Nel-| 
|son| ( |2006| and find that corrugation of the disc does not develop. 
Although these MHD simulations adopt a significantly lower res- 
olution than the pure hydrodynamic runs we have presented here, 
we note that hydrodynamic runs performed at low resolution still 
show the development of the instability even when the short ra- 
dial wavelength perturbations of the initial growth phase are not 
resolved. Instead, we find that the disc displays longer wavelength 
breathing and corrugation modes that become unstable and cause 
the disc to oscillate vertically in quite a violent manner. In mag- 
netised global disc models with dead zones whose vertical height 
covers ^2.5 scale heights, which support Reynolds stresses in the 
dead zone with an eff'ective value of a ^ 10~^, the development 
of these corrugation oscillations is observed in models that adopt a 
locally isothermal equation of state with q = -I. 

5.5 Thermal relaxation in models with T(R) 

We now consider the evolution of models where we relax the lo- 
cally isothermal assumption associated with the response of the 
fluid to perturbations. We evolve the energy equation in ([T]), and 
introduce thermal relaxation by integrating eqn. |T8] ). We adopt 
the equation of state P = (y - l)e, and set y = 1.4. The gas is 
assumed to be inviscid. Power-law profiles for the initial tempera- 
tures, T(R), and midplane density, Pmid(^), are adopted with q = -I 
and p = -1.5 in eqns. ^ and ([3]). The aim of these models is to 
examine the robustness of the vertical shear instability as a func- 
tion of the thermal relaxation time, Treiax, defined in eqn. (T8] l and 
expressed as a fixed multiple or fraction of the local orbital period. 
These runs are labelled T5R-0.01 - T9R-oo in Table 1. 

The evolution of the normalised perturbed kinetic energies for 
a number of models with relaxation times in the range < Treiax ^ 
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Figure 10. Time evolution of the sum of the (normalised) perturbed radial 
and meridional kinetic energy in discs where the temperature was initially 
constant on cylinders, as a function of the thermal relaxation time. Note that 
only the TReiax = and 0.01 cases show growth. 

oo are plotted in Fig. [To| It is immediately obvious that instability 
only occurs in either the locally isothermal case (Treiax = 0) and 
when Treiax = 0.01 orbits. All other simulations result in the per- 
turbed kinetic energy contained in the initial seed noise decaying 
with time. We note that the case of Treiax = ^ is directly compara- 
ble to a previous study on the adiabatic evolution of a stratified disc 
by lRiidiger et al.| ( [2002| ). The authors considered the hydrodynamic 
stability under the Solberg-H0iland criterion and also find stability 
in this case. 

Our results indicate that the vertical shear instability requires 
that the initial temperature profile of the fluid is re-established 
rather rapidly during dynamical evolution, at least for the equilib- 
rium temperature and density profiles adopted in these particular 
models. 

The requirement for near-isothermal evolution suggests that 
the vertical shear instability is most likely to operate in the opti- 
cally thin regions of astrophysical discs whose global properties 
are similar to those considered here. For example, the outer regions 
of protoplanetary discs lying beyond ~ 50-100 AU may be prone 
to this instability, provided that MHD turbulence is present at low 
enough levels that the instability is not damped by the turbulent 
viscosity. This seems to be a likely prospect given that low density 
regions may be stabilised by ambipolar diflfusion (Armitage'2011V 
It should also be noted, however, that the simple thermal relaxation 
model we employ does not capture the fact that the thermal evolu- 
tion time of a mode with radial wavelength Ar scales as ~ A^/D 
(where D is the thermal diflfusion coefficient), such that very short 
wavelength modes may remain unstable in optically thick discs. A 
reduction in spatial scales on which the instability operates, how- 
ever, will presumably affect the resulting turbulent flow and reduce 
the associated Reynolds stress and transport coefficients in the non- 
linear saturated state. 



5.6 Thermal relaxation in models with Ks(R) 

We now consider models in which the initial entropy function, 
Ks(R), follows a strict power-law function of radius given by 
eqn. ([6]), and the midplane density, Pmid(^), follows a radial power- 
law given by eqn. ([3]). For all models except one we adopt the val- 
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Figure 11. Time evolution of the sum of the normahsed radial and merid- 
ional kinetic energies in discs where the entropy function, KsiR), was ini- 
tially constant on cylinders. 



ues ^ = -1 and /? = 0, leading to the midplane Mach number, 
Almid. being constant at all radii. Our normalisation of Ks(Ro) sets 
Mmid = 20. We also consider a single model with p = -1.5 and 
5 = 0, so that there is no initial radial entropy gradient. The entropy 
in this case is normalised so that Almid = 20 Sii R = Rq = \. These 
runs are Hsted in Table 1 as KlR-0 to KlOR-0.01. 

We impose reflecting conditions at the meridional boundaries, 
and consider inviscid evolution. As described in Sect. |2.1.1| these 
models are convenient to implement numerically because analytic 
solutions can be obtained for the equilibrium density and velocity 
field through eqns. ^ and |T6| . As such, these models allow us 
to explore the vertical shear instability as a function of the ther- 
mal relaxation time, Tj-eiax. in discs where the initial distribution of 
temperature no longer follows a power-law function of cylindrical 
radius, but instead varies with both R and Z. We note that eqn. |T6| 
also demonstrates that a radial power-law in Ks(R) with 5 = -1 
gives rise to an equilibrium that varies with height Z at each ra- 
dius R, but adopting ^ = implies that no vertical gradient in 
exists. 

The normalised sum of the radial and meridional kinetic ener- 
gies is plotted in Fig. [TT] for thermal relaxation times in the range 
< Treiax ^ ^ local orbits. The model with 5 = and p = -1.5 
employed a relaxation time Treiax = 0-01 orbits, and is labelled as 

Treiax = 0.01, 5 = 0'. 

Interestingly, for the model with 5 = -1 and p = 0, all values 
of Treiax result in growth of the perturbed meridional and radial en- 
ergies, and only the strictly isentropic simulation with Treiax = ^ 
shows eventual decay over long time scales of ~ 100 orbits. It 
is noteworthy that this isentropic model is unstable according to 
one of the Solberg-H0iland criteria discussed in Sect.|3] and this is 
probably the cause of the initial growth in perturbed energy. The 
subsequent decay may arise because the instability causes the equi- 
librium disc to move to a stable state. Inspection of the perturbed 
velocity profiles in contour plots similar to Figs. |2]and[3] (not shown 
here) for all of the runs with p = and 5 = -1 shows that the previ- 
ously discussed characteristic perturbations with \kR/kz\ » 1 grow 
initially, even for the isentropic disc model with no thermal relax- 
ation. All models with finite thermal relaxtion time show long term 
growth in perturbed energy, although in the case of Treiax = 10 or- 
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Figure 12. Time evolution of the perturbed meridional + radial kinetic en- 
ergy (normalised by the total energy in keplerian motion) for the full 3D 
simulation T1R-0-3D. 

bits the growth time is very long indeed. The instability displayed 
by the remaining models with Treiax ^10 orbits shows that the re- 
quirement of very rapid thermal relaxation observed in the models 
presented in Sect. |5. 5 [ actually depends on the detailed temperature 
and density structure of the disc. It is clear that disc models exist for 
which thermal relaxation times in the range < Treiax < 10 orbits 
lead to the growth of the instability, and as such its range of applica- 
bility in the study of the dynamics of astrophysical discs is probably 
broader than suggested by the results presented in Sect. |5.5| A full 
exploration and understanding of the range of applicability, how- 
ever, is beyond this paper, and will require a dedicated and detailed 
future study that accounts for the thermal evolution of the disc with 
greater realism. 

Turning to the run with p = -1.5 and 5 = 0, we see that 
the initial perturbation energy decreases for the first ~ 200 orbits 
before increasing again. The expectation is that this model will not 
display the vertical shear instability, and inspection of velocity con- 
tour plots (not shown here) confirms that the characteristic pertur- 
bations with \kR/kz\ » 1 do not appear in this case. These velocity 
plots, however, indicate that over secular time scales sound waves 
are generated close to the meridional boundaries, and this appears 
to be the reason for the up-turn in the perturbed kinetic energies 
seen in Fig. pT] after 200 orbits. 

5.7 Non-axisymmetric model 

We now consider briefly the evolution of a non-axisymmetric 
model T1R-0-3D, in which the azimuthal domain covered 7t/4 
radians. This model is the 3D equivalent of model TlR-0. The 
simulation was performed using the nirvana code, and the model 
was set up using equations ([3]), ^ and |T3] l, with values p = -1.5 
and q = -1. The velocity field was seeded with noise (amplitude 
0.01c J. Details of the model are given in Table[T] 

The total normalised kinetic energy (meridional + radial) ver- 
sus time is displayed in Fig. [12] Comparing with the equivalent 
axisymmetric plot (Fig[6| we see that the evolution is similar, with 
growth in the perturbed energy occurring after ~ 10 orbits, and sat- 
uration at a value of a few xlO~^ beginning to occur after ~ 100 
orbits. 

In Fig. [13] we plot the time evolution of the volume averaged 
Reynolds stress normalised by the mean pressure in the disc. This 
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Figure 13. Time evolution of the volume averaged Reynolds stress (nor- 
malised by the mean pressure) for the full 3D simulations T1R-0-3D 



quantity is computed as follows. We define an azimuthally aver- 
aged Reynolds stress T^{r, 6) obtained by averaging the quantity 
p6v^6vtp over azimuth. Here 6v^ and dv^p are the local radial and 
azimuthal velocity fluctuations. We also define a density- weighted 
mean pressure as a function of r, P(r), obtained by averging over 
6 and 0. We define a local value of the Shakura-Sunyaev stress pa- 
rameter a{r,6) - T^(r,6)/P(r). The simple average of a{r,6) over 
r and 6 is the quantity plotted in Fig. [13] 

Although rather noisy, we see that the normalised stress ap- 
proaches average values ~ 6 x 10""^ by the end of the simulation 
(and appears to be still growing at this point). The spatial distribu- 
tion of Qf(r, 6), time averaged during the last 10 orbits of the run, is 
shown in Fig. [14] Here we see that local values of the stress reach 
~ 2 X 10~^, indicating that the vertical shear instability generates a 
quasi-turbulent flow capable of supporting significant outward an- 
gular momentum transport in astrophysical discs, given favourable 
conditions for its development. 

The upper panels of Fig. [15] show contours of the perturbed 
density, 6plpQ in a slice parallel the meridional plane at three diff'er- 
ent times during the simulation, showing similar features to those 
presented for the 2D-axisymmetric simulation in Fig. [5] Perhaps 
more interesting are the lower panels of Fig. [15] which show the ac- 
tual density p in the {R, 0) plane located at the disc midplane. Here 
the development of spiral density waves may be observed, similar 
in morphology to those that arise in discs where turbulence is driven 
by the MRI (e.g. Papal oizou & Nelso n 2003 ). The 3D simulation 
presented here suggests that if the appropriate conditions prevail in 
astrophysical discs, the vertical shear instability may lead to a tur- 
bulent flow capable of supporting significant angular momentum 
transport. 



6 THEORETICAL CONSIDERATIONS 

The GSF instability appears by rendering the inertial modes of a 
rotating atmosphere unstable. The original analysis in Goldreich 
& Schubert (1967, GS67 hereafter) demonstrated the possibility of 
this instability by performing a point analysis at a given location in 
a stellar radiative zone away from the equator, equivalent to consid- 
ering a location away from the midplane in a disc. In this section 
we examine the mathematical structure of the instability by further 
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Figure 14. Spatial distribution of the time and horizontally averaged 
Reynolds stress (normalised by the mean pressure at each radius) for the 
model T1R-0-3D. 



extending previous analyses, including those of Urpin (2003) and 
Arlt & Urpin (2004), by relaxing the point assumption. 

We shall focus on disturbances which are locally isothermal. 
For the sake of completion of this important discussion we redo 
the original point analysis of OS 67 in Appendix A, but without in- 
troducing the Boussineq approximation. Denoting cr as the growth 
rate we find that the inertial mode response is roughly given by 



-Kl{clkl+Nl) + 2^1^clkrKf 

clikl+kD + Kl+Nl 



(26) 



in which cq is the reference sound speed, Kq is the epicyclic fre- 
quency which, for a keplerian disc, is given by Qq. the local kep- 
lerian rotation rate at the point in question, and ky are the corre- 
sponding vertical and radial disturbance wavenumbers respectively. 
The local Brunt- Vaisaila frequency is A^o and dV/dz is the vertical 
gradient of the mean keplerian flow. This quantity typically scales 
on the order of magnitude of (q/2)Qo(Ho/Ro) where q is the same 
exponent of the radially varying isothermal sound speed discussed 
in Section 2. Supposing for this discussion that A^o is negligible it 
follows from this expression that if Hq/Rq is small then instability 
can only happen if the radial wavenumber conspires to be corre- 
spondingly large. In that limit the above expression implies 



2^0 



LdV 



Kdz 



(27) 



indicating that instability is possible if k^/K ~ 0(qHo/Ro). The 
analysis of Arlt & Urpin (2004), for example, also similarly in- 
dicates that for the same rough conditions the growth rate ought 
to scale as 0(qQoHo/Ro). The simulations we have performed are 
consistent with this tendency where the radial length scales of the 
emerging structures are significantly shorter than the vertical ones 
with growth rates of the instability on the order of 4 orbit times for 
Ho/Ro ~ 1/20 and^ = -1. 

Our goal is to develop a better physical understanding of the 
processes responsible for this instability beyond invoking Solberg- 
H0iland criteria. In this respect we notice from Figure 1 that radial 
velocity fluctuations are considerably smaller in magnitude than the 
corresponding meridional velocities. 

With these clues in mind, we show in the following how 
the processes involved in bringing about the instability is largely 
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Figure 15. Perturbed density, Sp/p in the merdional plane at = ;r/8 (upper three panels) for the 3D simulation T1R-0-3D. Note that we have effectively 
stretched the grey-scale by plotting the quantity sign((5p) x |(5p/pp/^ in the upper panels. The lower panels show the relative density perturbations 5p/p at the 
disc midplane. No grey-scale stretching has been applied to these lower panels. 



anelastic and radially geostrophic - by the latter expression we 
mean to indicate dynamics which are in constant radial force bal- 
ance between Coriolis effects and pressure gradients. Furthermore, 
despite the varied simplifications we make to expose the essence of 
the physical process, the fundamental equations describing the re- 
sulting linearised response remain inseparable in the radial and ver- 
tical coordinates. This means that the only recourse in establishing 
any insight is through a further approximate solution of the result- 
ing reduced equations. We find that the solution indicates that for 
given parameters describing disturbances the instability appears in 
pairs, as opposed to appearing individually as indicated by ( [27| . 
Although we have not proved it in this study, we conjecture the 
powerful driving of the instability in the simulations may be, in 
part, caused by this feature. 



6.1 Equations of motion revisited and steady states rederived 

The equations of motion for axisymmetric inviscid dynamics in a 
cylindrical geometry are given by 



d_ 

dt 

d 



u- 



d_ 
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w- 



dz 



~ R 

uv 
w 



1 dc^p 
~p~dR 
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1 dc^p 

~p~dz' 



do 
dR' 



do 

dz' 



(28) 



(29) 



(30) 



Note that the (R, 0, Z) velocity components are given here by 
(U, V, W). We dispense with the subscripted scheme (v^, v^, vz) used 
in previous sections in order to simplify the notation. The corre- 
sponding equation of mass continuity is 

dp 1 dRpU 



dt ^ R dr 



dpW 

~dz 



As mentioned above, we focus here on dynamics that are locally 
isothermal with an infinitely short cooling time (Treiax 0)- This 
then is to be considered in the context of simulations TlR-0 to 
T4R-0 summarized in Table 1. Reciting therefore from Section 2: it 
means that the square of the sound speed is given by c] - cliR/RoY 
where Rq is the fiducial reference disc position and cq is the scaled 
sound speed at that point. The gravitational potential emanating 
from the central object is O = -GM/(R^ + Z^f'^. 

The general equilibrium state solutions are found in eqns. 
(12)-(13) but, as we mentioned earlier, perturbations superposed 
on this base state are difficult to analyse because the resulting equa- 
tions are fundamentally inseparable so that a typical normal-mode 
analysis is out of the question. In order to facilitate some kind of 
tractable analysis we make the one and only approximation here: 
the radial and vertical gradients of the potential O are expressed in 
terms of their corresponding first order Taylor Series expansions, 
i.e. 



do 
dR 
do 

dz 



GM 
GM 



(If 



= 0. 



(31) 



in which Qq = (GM/R^y^^ is the reference keplerian rotation rate 
at radius Rq. The midplane density is chosen to be of the form 
Pmid = PoiR/RoY where po is the reference density and where p is 
an arbitrary index (as referenced earlier). In the following analysis 
it will be convenient for our discussion to refer to the natural loga- 
rithm of the density instead of directly to the density itself, thus, we 
define 11 = Inp. Because we shall be concerned with perturbations 
around the steady states implied by the above equations we shall 
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represent these states by overbars. As such, we have that 

. 1 Z2 / 1 

n^lnp.,--- - 



1 + s{R) - 



qlHo 
A\Ro 



(32) 



,(33) 



where Hq = cq/Qq is the local vertical scale height as referenced 
near the end of Section 2. The non-dimensionalised scale height 
is accordingly given via the relationship TY^ = (R/Rq)^'^^ also as 
found in Sect. 2. The nondimensional quantity s is given by 



[Roj \Ro 



^0 TTp. 1" ^0 " 



= (p + q) 
where = (R/Roy. 



dR 

l+q 



dR 



6.2 Linearised perturbations and non-dimensionalisation 

We introduce perturbations by writing for each dependent quantity 

and inserting these into the governing equations ([28]l-([37]). Lin- 
earising results in the expressions 



dv ,idRV ,dv 

^^''R^^'^dZ 
dw' 

~dt 



0, 



, dIV 
'~dR 



y dW 

'~dz 



(34) 



and 



dU 



dU 



(35) 



dW _ 1 dRu' dw' 

"df ~ ~R~dR' ^ ~dz ~^ dR~^ dz' 

It will now be made more transparent if we non-dimensionalise the 
above equations according to the quantities appearing. We see that 
a natural time unit is given by the keplerian rotation time The 
radial and vertical length scales are naturally scaled by Rq and Hq 
respectively. Thus we write for the these quantities 



R Ror, 



(36) 



where t, r, z are the corresponding non-dimensionalisations of the 
independent variables representing respectively time, radius and 
height. It is very important to note that the radial and vertical length 
scales are disparate with respect to each other by a factor of Hq/Rq. 
Because in all of our simulations this ratio is quite small (~ 0.05) 
we shall treat this ratio as one of our "small parameters" and for- 
mally represent it by 6 = Hq/Rq (not to be confused with £(R) 
defined earlier). This disparity must be kept in mind when the scal- 
ings invoked to recover the GSF instaiblity are formally made in 
the next section. 

Judging from the dynamics observed in the simulations the 
structures appearing tend to be radially and vertically constrained. 
These spatial constraints (especially the radial confinement) in- 
dicate that perturbation velocities ought not to exceed the sound 
speed (at least initially). This is typical of the scalings frequently 
used to derive equations appropriate to the dynamics in a small 
box of a disk ( [Goldreich & Lynden-Bell| 1965 [[Umurhan & Regev| 
[2004[ ) although we are not, technically, considering the dynamics 



on such small scales yet. In sum, therefore, we scale the dependent 
perturbation velocities by 



CqU, 



CqV, 



W CqW, 



where, as before, u,v,w represent the corresponding non- 
dimensionalised component velocities in the radial, azimuthal and 
vertical directions. Therefore, the perturbation equations now take 
on the following more transparent appearance 
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(37) 
(38) 
(39) 



— - -ei- + —\-6—u- — -w— (40) 
dt \r dr I dr dz dz ' 

where the non-dimensionalisation of the mean azimuthal flow V is 
given in terms of the other redefined variables 

1 



2'H2' 



V 



QoRo 



-3/2 



(l+£(r)+^6Vr- 



(41) 



with s(r) = 6^(p + qy^^"^ and 'K^(r) = r^+^. We have written ( 38 1 in 



a seemingly curious way: the last term on the RHS of that equation 
is the product of the vertical gradient of the mean azimuthal flow 
term, i.e. -wdvjdz. We have chosen to write it out explicitly in 
order to bring to the fore the leading order scaling that sits in front 
of it as it will efl'ect how we proceed toward the reduced model (see 
the next section). 



6.3 Asymptotic scalings and resulting reduced equations 

The linearised equations of motion ([37])-([40| are, despite our efl'orts 
to simplify, still inseparable between the r and z variables. In order 
to proceed asymptotically we must make further scaling choices. 
These are guided by both the results of the numerical solutions as 
well as by the discussion at the beginning of Section 6. At this stage 
we shall list them once again: 

(i) As estabilished by GS67, Urpin (2003) and Arlt & Urpin 
(2004), growth rates scale as ~ ^Qq^o/^o which is, in our non- 
dimensionalised time units, ~ eq, 

(ii) For growth rates slow on the timescale of the local disc rota- 
tion, the emerging structures have radial dimensions (4) consider- 
ably smaller than the corresponding vertical dimensions (4)- That 
is to say, for eq <^ I the scaling analysis of both GSF67 and Urpin 
(2003) indicate that 4/4 ^ 1 

(iii) The numerical solutions also clearly indicate that during the 
growth of the instability the radial velocity fluctuations are signifi- 
cantly smaller than the corresponding meridional velocity fluctua- 
tions (see Figure 1). 

In the following we describe scalings of ([37]l-([40| that sim- 
plify them into a set which is both more transparent and more 
amenable to further analysis while retaining the essential physi- 
cal processes involved in the instability. We assume that 6^1 
and treat q as an order 1 quantity (although a more general analy- 
sis can be done without this a priori assumption achieving, in the 
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end, much of the same results discussed hereafter). Furthermore we 
consider the analysis around the fiducial radius r = 1. Since interest 
is in radial scales that are much smaller than the vertical scales and 
recalling that the r scales dimensionally represent physical length 
scales that are longer than the dimensional vertical scales z by a 
factor of 6"^= Ro/Ho ~ 20) we consider radial disturbances 

r - 1 = e'^x, 

where x is order 1 . We leave the z scales untouched as these are the 
de-facto reference scales of the analysis. Because the growth rates 
are long by a factor of we introduce a new long-time variable r 
given by 

t = Te-\ 

With this long-time scale assumed we find that in order to bring 
about non-trivial pressure balancing with the inertial term in the 
vertical momentum equation it must follow that the pressure fluc- 
tuations must relatively scale by € as well. This can be easily sur- 
mised by examining ([39} and noticing that for w order 1 and the 
time derivative scaling as order 6, that the only way balance occurs 
is if the pressure is correspondingly small by a factor of 6. This 
means introducing a new pressure fluctuation reflecting this scal- 
ing through 

where ft is the scaled pressure. 

Finally as we have just intimated, in addition to the vertical 
velocity being order 1 we assume that the azimuthal velocity fluc- 
tuations are also unsealed (i.e. remaining order 1) in accordance 
with our numerical observations. We note here that scaling v to be 
order 1 is also consistent with the pressure scalings assumed be- 
cause it leads to a balance between the radial pressure gradient and 
the Coriolis term in ([37]). 

However, we suppose that the radial velocities are small in 
comparison to the other velocity components and we propose that 
its relative smallness is similar to the pressure field's scaling, i.e. 

u - €U, 

where u is the correspondingly scaled radial velocity. Applying 
these scalings assumptions to eqns. ([37]l-([40| results in the follow- 
ing equations at lowest order, 



(42) 

qzw, (43) 

(44) 

- zw, (45) 



with corrections to the above equations appearing at order 6^. In 
this form these reduced equations contain insight with regards to 
two very important physical implications. The first of these follows 
from the interpretation of eqn. ([42} which says that the dynamics 
of the instability occur under radially geostrophic conditions, that 
is to say, that the processes develop under conditions in which ra- 
dial Coriolis eff'ects balance radial pressure gradients. The second 
observation is that the linear dynamics are anelastic rather than in- 
compressible in character. By this we mean to say the following: 
since on these radial/vertical length scales the mean (scaled) den- 
sity profile has the form p = e~^ eqn. {45} may be equivalently 



written as 
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dx dz 



(46) 



= 


2v- 


du 

dx ' 


dv 


1 

"2^ 


1 

"2 


dw 


dii 










= 


du 


dw 


dx 


'di 



The fact that the dynamics here are not incompressible in the usual 
sense is perhaps less surprising given that vertical stratification is 
non-negligible under these spatial constraints - had we been inter- 
ested in vertical scales that were equally as short as the radial scales 
then stratification would not figure prominently. 

The other two equations describing the vertical and azimuthal 
momentum balances retain their inertial terms and are largely un- 
aff'ected (directly) by these scalings. 

Before analyzing the solutions to these equations it is impor- 
tant to keep in mind that the essential eff'ect giving rise to the in- 
stability is present in the guise of the final term on the RHS of 
eqn. ( [43} . Additionally, in reflecting upon these equations, it should 
be kept in mind that u and ft indicate real quantities that are intrin- 
sically smaller (but not zero) compared to the other terms. 



6.4 Approximate solutions and double instability 



The reduced system ( |42|45| l may be combined into a single equa- 
tion for the pressure perturbation ft: 



_^#ft__#ft^^ m 

dr^ dx^ dz} \ ^ dxj^ dz 



(47) 



We note here that a point analysis of ( [47} , i.e., assuming z = Zo is 
fixed and making a wave ansatz and proceeding similarly to Ap- 
pendix A recovers the content of the asymptotic growth rates con- 
tained in eqn. ( [27} , indicating the consistency of the scaling argu- 
ments we have exploited to get to this point. 

Nonetheless, examination of this equation, although appearing 
quite simple in many respects, shows that it is also fundamentally 
inseparable owing to the qd/dx term on the RHS of the expression. 
A more concerted future analysis is necessary to develop proper 
solutions of these equations subject to proper boundary conditions 
on both the vertical and radial boundaries of the system. However, 
we may proceed analytically by the following rationale. 

Keeping in mind that that equation ( [32} says that the basic 
state density profile p has the form ~ e~^ we can develop solu- 
tions of these equations in which we require that 

(i) The pressure fluctuations ft are zero on the radial boundaries 
located at x = +A and, 

(ii) The kinetic energy in the fluctuations decay as z +oo. 
Since the kinetic energy involves terms that are appear as pif,pv^ 
and pw^ then this boundary condition is satisfied even if the velocity 
fields show algebraic growth as z ^ +oo. 

The first boundary condition is not the same as the no-normal flow 
boundary conditions of the simulations but we have checked and 
found that the resulting growth rates and qualitative results are un- 
changed when implementing those conditions instead (see footnote 
later). This second condition is not outright unphysical but it, nev- 
ertheless, does not mirror the conditions present in the simulations 
and we must keep this in mind when we analyse the results later 
on. We assume all fields have normal-mode form in time 

ft(x,z,0 = tL(x, z)e'^ + C.C. 

where s is the growth rate. This form is also assumed for the other 
variables u, v, w. As a side note we point out that with this solu- 
tion ansatz inserted into eqn. ( 44 1 it follows that sw = -^^ft- The 
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equation for fl becomes 
2#n _ 



d \ dii 



dx dz 



(48) 



We proceed to develop solutions of ( |48) ) guided by the methods 
used to develop solutions to Hermite and Gegenbauer differential 
equations (|Morse & Feshbach 1953 , Lin 2012 ). With m a positive 
integer we assume a tractable non- separable solution of the follow- 
ing form 



U(x,z) = J]UjM)^, 



(49) 



where the radial functions TLjmix) are yet to be determined. [^Insert- 
ing this solution ansatz into ( [48] l reveals that for a given m and for 
each index j - ■ ■ ■ ,m - 4,m - 2 the following equations must be 
satisfied 



dx^ 



dx 



(50) 



together with the corresponding "top" equation, i.e. the ordinary 
diff'erential equation for j = m, 



, + qm- 



dx 



(51) 



The index j starts from (1) depending if m is even (odd). In re- 
lation to the modes observed in our numerical solutions: breathing 
modes correspond to even values of m while corrugation modes 
correspond to odd values of m. We also note that each function 
tijfn{x) must individually satisfy the boundary condition IIy,„(+A) = 
for each j. This is a rather cumbersome task and it is not our in- 
tention here to fully develop detailed solutions. Rather, we seek to 
determine the value of the growth rate s which may be calculated 
by solving the top equation ( [5T] ) subject to the boundary condition 
n^OT(+A) = 0. Elementary analysis of this equation indicates solu- 
tions are of the form 



COS kx. 



where n - 1,2, 



(52) 



where is an arbitrary constant. The requirement for the consis- 
tency of this solution follows from inserting \52) directly into 
which, after factoring out limm reveals the growth-rate relationship 
for s 



s + 



m 



2 2 



which has solution 



m 



(-1 + Vi - q^k^) ■ 



(53) 



(54) 



It follows from an examination of the above equation that for 
< < 1 the two roots possible for are both negative which 
implies that the four modes are all oscillatory. However, instability 
is possible once qk > \ since it implies that is now a complex 
number with a non-zero real part. Moreover, because of the sign 
convention, for qk > I there are always a pair of unstable modes 
and a pair of stable modes. In other words we always have in that 
case the four possible combinations: s = ±(\Sr\ ± i\si\) for some real 

^ For example, when developing solutions to the Hermite differential equa- 
tion of some integer order, the constant coefficients in front of each power 
of z are are determined through a standard recursive procedure in descend- 
ing powers of z. The same is done here except that the recursive procedure 
generates coupled ordinary differential equations in x. 



values of Sr, Si > 0. This quality holds also for no-normal boundary 
conditions at x = +A.[^For unstable values (\kq\ > 1) the expres- 
sions \Sr\ and \si\ are given by 



\Si\ = 



^J\kq\ + l. (55) 



\k\ ' \k\ 

An elementary examination of the growth rate Sr shows that the 
wavenumber corresponding to maximal growth occurs for \kq\ = 2. 
If we let /Imax denote the corresponding wavelength in dimensional 
units then this wavelength of maximal growth is 

, - ,? ("o^ 

where, for our expreriments corresponding to ^ = -\,HqIRq - 
0.05, implies radial scales on the order of /^max ~ 0.008/?o- These 
scales are approximately the length scales observed at the early 
stages of growth in both Figures 2 and 3, especially near the Rq, = I 
boundary. Similarly, the growth rate of the fastest growing mode 
denoted by 9i^(5max) = \q\ ^mllll and, in dimensional units, this 
implies a maximum growth rate CTj^ax 

o-max = i=r — Qq = ym\q\—- — orbit . 

2V2 \^o/ V2\^o/ 

For the simulations TlR-0 to T4R-0 (Table 1) this means a maxi- 
mum growth rate given by (Tmax ~ 0.11|^| ^/m orbit" ^ It should be 
kept in mind that these growth rates are for the individual velocity 
fields. However, the growth rates reported in Sect. 5 are for the per- 
turbation kinetic energies which depend upon the square of the ve- 
locity fields. As such, the linear theory therefore predicts the max- 
imum growth rates /6>r the kinetic energy perturbations to ~ 2o-^ax 
which, if applied to our numerical experiments, would be 

~ 0.22M orbit"^ 

This is significant since we reported there that the growth times 
were measured to be ~ 0.24 orbit"^ The theory predicts: for the 
fundamental corrugation mode m = 1 a growth rate of ~ 0.22 
orbit" \ while for the first breathing mode m = 2 a growth rate 
of ~ 0.31 orbit" ^ This suggests that the growth rates observed in 
the early stages of growth of simulations TlR-0 to T4R-0 is a con- 
volution of these two fundamental modes of the system. 

The theory, as it stands, does not select for a preferred value 
of m as the growth rates all scale as ^/m suggesting curiously that 
for larger values of m one should expect faster growth rates. This 
appears to be in conflict with the original point analysis of GS67 
and others. In particular, since m roughly indicates the number of 
nodes in the vertical direction one can treat m as a proxy for the 
vertical wavenumber k^. But according to the criterion obtained via 
the point analysis ([27]), increasing the vertical wavenumber ought 
to promote stability. This conflict of interpretation needs to be re- 
solved in follow up work, however, we do offer two hypotheses 
regarding this matter: 

(i) Since the instability appears to be global in character (as 
evinced by the numerical solutions) that the classical point analy- 
ses serve only as a useful guide toward indicating the possibility of 



^ In this case, after inserting the solution ansatz into the top equation p3) 
indicates that the no-normal boundary condition is enforced by requiring 
s^dxti-mm - qmtimm = at X = ±A. What emerges is a more complicated 
relationship between s and k, m but resulting in the same qualitative features 
discussed for fixed pressure conditions. A more thorough exposition awaits 
a subsequent study. 
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instability only. Additionally, since the reduced equations derived 
here are more global in reach, it would be logically more appropri- 
ate to be guided by their content in better describing the emergence 
of the instability and/or, 

(ii) The implementation of the decaying kinetic energy bound- 
ary condition for z ^ ±00 may be somehow responsible for the 
feature regarding m and the mode's stability character. It could be 
that high m modes may become stable if one, instead, imposes no 
normal flow conditions at some finite height away from the disc 
midplane. Of course, to do this requires non-trivial solution forms 
that are far more complicated than the solution ansatz assumed in 



7 DISCUSSION AND CONCLUSION 

We have presented the results from a large suite of 2D- 
axisymmetric and 3D simulations, using two independent hydrody- 
namic codes, that examine the stability and dynamics of accretion 
disc models in which the temperature or entropy are strict functions 
of cylindrical radius. Such thermal profiles lead to equilibrium an- 
gular velocity profiles that depend on radius and height. We find 
that these disc models are unstable to the growth of modes with 
\kR/kz\ » 1 when the thermal evolution time is comparable to or 
shorter than the dynamical time. This mode pattern is expected for 
systems undergoing the Goldreich-Schubert-Fricke instability. The 
potential for this vertical shear instability to apply to accretion discs 
has been investigated previously by|Urpin & Brandenburg |(T998] ), 
|Urpin| | |2Q03l l and |Arlt & Urpin| ( |2004| ) using the Boussinesq ap- 
proximation, and we have confirmed this with our own analysis that 
applies to fully compressible flows. The nonlinear simulations in- 
dicate that evolution of the instability leads to the dominant radial 
wavelength increasing from a fraction of the scale height during 
early times to being comparable to the scale height at later times. 
Instability also appears to generate disc vertical motions that corre- 
spond to vertical breathing modes at early times, but at later times 
the dominant motion corresponds to corrugation of the disc caus- 
ing the disc midplane to oscillate around its initial location at all 
radii (at least in axisymmetric flows). A full 3D simulation com- 
puted using a locally isothermal equation of state indicates that the 
instability generates a turbulent flow in its nonlinear saturated state, 
leading to non-negligible transport of angular momentum through a 
Reynolds stress with corresponding viscous alpha value of ~ 10~^. 

The requirement for quite rapid thermal evolution of the disc 
for strong instability to ensue suggests that the outer regions of pro- 
toplanetary discs may be places where this instability operates be- 
cause of the low density and optical depth there. We caution that 
the thermal evolution time required for instability may depend on 
thermal gradients, such that instability can operate for longer cool- 
ing times in discs with steeper temperature/entropy profiles, but it 
seems unlikely that instability can be sustained if the local thermal 
time scale greatly exceeds a few orbital periods. We have shown 
that the instability is damped in viscous flows with dimensionless 
kinematic viscosity v ^ 10~^, so these outer disc regions would 
need to be stable against the MRI for the instability to operate (the 
instability is not observed in global simulations of discs support- 
ing fully developed MHD turbulence). It is likely that the MRI is 
quenched by ambipolar diffiision in the outermost regions of pro- 
tostellar discs ( |Armitage|2011| ), so these may provide low density, 
magnetically inactive regions where the vertical shear instability 
can operate efficiently. 

The vertical shear instability is an example of a baroclinic in- 



stability because the required vertical shear in a disc only occurs if 
the pressure is a function of both density and pressure, P(p, T). An 
important question that we have not addressed in this work is how 
a disc would evolve that is subject to both the vertical shear insta- 
bility and the subcritical baroclinic instability (SBI) studied by |Pe-| 
Iter sen et a'L] (|2007^ and'L esur & Papaloizou| ( [20 10). The conditions 
required for the SBI to operate and be sustained strongly are a radial 
entropy gradient and a fairly rapid thermal relaxation time, and we 
have shown that the vertical shear instability operates under these 
conditions also. Both |Petersen et al.| ( |2007| ) and |Lesur & Papaloizouj 
P010| ) report that a strongly sustained instability is obtained for 
thermal time scales close to the local orbital period, leading to a 
highly complex flow in which long-lived vortices are formed, and 
an effective viscous stress of a ~ few xlO~^ is maintained by the 
Reynolds stress in compressible flows. The 3D simulation we pre- 
sented in Sect. |5.7| used a locally isothermal equation of state, and 
so is not subject to the SBI, but it seems likely that the combined ac- 
tion of the two instabilities in a disc with longer thermal evolution 
time and appropriate entropy gradient will generate a complex flow 
containing long-lived vortices and vertical motions that correspond 
to corrugation of the disc, accompanied by a Reynolds stress that 
leads to efficient outward angular momentum transfer. It is worth 
noting that the vertical shear instability is linear, whereas the SBI 
is a finite- amplitude instability, so it is possible that the SBI may be 
stimulated by perturbations generted during the development of the 
vertical shear instability. We will present a study of these two insta- 
bilities operating in tandem in a future publication to explore these 
hypotheses. Given the role that vortices may play in the trapping 
of solids during planet formation (e.g. [Barge & Sommerial|1995[ 
[Klahr & Bodenheimer|2003l ), this is clearly an important issue for 
further investigation. 
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APPENDIX A: RECITAL OF GS67 WITH LOCAL 
ISOTHERMAL ASSUMPTION 

Owing to the importance of the original study presented by GS67, 
we repeat here the calculation contained in that work wherein we 
assume the equation of state of the gas is locally isothermal with a 
temperature profile varying in the nominal radial direction. We also 
extend that analysis to include the eff'ects of compressibility and we 
do not assume outright that the hydrodynamic flow is incompress- 
ible. 

The original analysis was performed in a local rotating refer- 
ence frame (the original intention being to examine the possibility 



of such instabilities in the interiors of rotating stars). The equations 
of motion considered in that work are akin to the local shearing 
sheet approximation familiar to accretion disc theory ( [Goldreich &| 
[Lynden-Bell 1965 ). Since this is our concern here we shall use that 
as our starting model set. The equations for axisymmetric dynamics 
in that environment are therefore 

du del 



dv 
dt 



+ IQqu 



0, 



dz 



(Al) 



dw 
dt 

dli _|_ _|_ Q 
dt dx dz 

where the expressions appearing are consistent with those devel- 
oped in the body of the text recalling especially that IT = In p and 
Qo = constant (see Sect. 6). The remaining expressions are de- 
fined again for convenience below. The only diff'erence is that we 
use lower case letters for the radial, aziumthal and vertical veloc- 
ities (i.e. u, V, w) in order to distinguish them from the ones used 
to describe dynamics in a cylindrical geometry examined in the 
main text of our study. Typically speaking, the vertical component 
of gravity is vertically varying and is given by g = Q^z - but /or this 
analysis it is treated as a constant. 

As done in GS67, one can do a point expansion (see the dis- 
cussion of GS67 right prior to eq. 17 of that work) around any nom- 
inal level z = Zq. We start by considering the mean states which we 
represent with overbars. Thus, 

-2QoV = -c?—--^ 







dx 



dx 



Note that the mean azimuthal flow state v has been decomposed 
into a keplerian piece (the term -(3/2)Qo-^) plus a deviation about 
that state V, i.e. v = -qQqx + V. As before, is the sound speed 
(implicitly a function of x because of the radial dependence of the 
vertically isothermal temperature profile). In steady state we find 
that 

dz 

and, most importantly, the mean gradient of the azimuthal flow is 
dV g dlncl 



(A2) 



dz 2Qo dx 

Perturbations of IT and all other variables around their reference 
means states are introduced with prime notation, e.g. IT ^ H + 11', 
etc. We assume an isothermal equation of state for the perturbations 
as well. Linearised perturbations of the equations of motion reveal 



du' 
~dt 



■ 2Qov' 



dV_ 
'dt 



+ ^lo[\l2 + V,)u' + V, 



dw' 

~dr 

dW ,(9n du' dw' 

-IT- +W — +W — + — + ^— 

dt dx dz dx dz 



y dW 

'~d^ 



, dW 

'~d^ 



0. 



(A3) 



Note we have utilised the shortened notation = dV/dx. To em- 
phasize: the analysis we carry out here pointedly departs from that 
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done in GS67 in two respects: (i) we assume an isothermal equation 
of state for the disturbances and (ii) we allow for compressibility 
(cf. eq. 29 of GS67). Since g is constant, the above equations are 
easily combined into a single one for 11' yielding, 



+ 



m 

dx 



d_ 

dx 



m d_ 

dz dz 



df 

: (A4) 



where the epicyclic frequency is represented by kq and is related to 
the steady state quantities by 

Paraphrasing directly from GS67 (right prior to eq. 17 of that 
work): The next step is to expand the unperturbed variables and 
their derivatives in Taylor series about some point x = xq and 
z = Zq. The latter restriction is not as arbitrary as it might seem 
- it just makes the calculation more transparent. Discarding terms 
of order (x - xq)Ixq and (z - zo)lzo the perturbation variables may 
be expanded in plane waves of the form ~ ^lioji+k^x+k^z) ^ revealing 
the dispersion relationship 

in which Cq is the sounds speed at the point in question and where 
Kq O^q. As is standard for problems of atmospheres, we can recast 
the vertical wave dependence to have a complex character (indicat- 
ing a basic wavey pattern) by defining ^ - ig/(2cl) which 
renders the dispersion relation into the form 



.[cl(kl+kl) + Kl+Nl]co'- 



■ 2QoV,clk,k, + Kliclkl + A^o^) 



0, 
(A5) 



This equation has imaginary solutions for oj if (i ) C < with B > 
or (ii) - 4C < 0. The GSF criterion ( |A6| ) is essentially the 
condition that C < 0. Condition (ii) is for acoustic-inertial modes 
to be unstable. This amounts to 

[clikl + kl) + 4+ N'of - 4 [-2QoV,clkA + ^liclkl + A^o')] < 0. 

(A7) 

A detailed examination of the behaviour of this condition shows 
that it does not occur for reasonable values of the parameters lead- 
ing us to conclude that the acoustic modes (per se) do not become 
unstable as well. 



where we have defined the Brunt- Vaisaila frequency 



Solutions for o) are easy to obtain and write down. It is more in- 
structive, however, to assess the stability characteristics straight 
from an analysis of \A5\ itself. 

The classical GSF instability is one in which the modes pass 
through zero frequency before becoming unstable and they de- 
scribe the influence that the vertical shear gradient has upon iner- 
tial modes. In our dispersion relation this amounts to the instability 
condition 



- 2^1oV,clk,k, + Kliclk] + Nl) < 0. 
When c^kl » the condition reduces to 



(A6) 



-2^qV^ + kIy 

which is essentially the first term appearing in eq. 33 of GS67 (mis- 
printed). Also, this condition is identical to that found in Urpin 
(2003), eq. 20. These latter correspondences follow from realizing 
that 



1 IdR^Q. 



2^^o ^ Rl 



dR 



R=Rq 



and, thus, recovering the classical Goldreich-Schubert condition. 
The dispersion relation \A5\ has the general form 



Boj^ + C = 



where B and C are obviously identified with the terms in \A5\ . 



